home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / srcuc.zip / FFT.C < prev    next >
C/C++ Source or Header  |  1991-12-20  |  60KB  |  1,744 lines
  1. /* -*-C-*-
  2.  
  3. $Header: /scheme/users/cph/src/microcode/RCS/fft.c,v 9.30 1991/12/20 22:49:07 cph Exp $
  4.  
  5. Copyright (c) 1987-91 Massachusetts Institute of Technology
  6.  
  7. This material was developed by the Scheme project at the Massachusetts
  8. Institute of Technology, Department of Electrical Engineering and
  9. Computer Science.  Permission to copy this software, to redistribute
  10. it, and to use it for any purpose is granted, subject to the following
  11. restrictions and understandings.
  12.  
  13. 1. Any copy made of this software must include this copyright notice
  14. in full.
  15.  
  16. 2. Users of this software agree to make their best efforts (a) to
  17. return to the MIT Scheme project any improvements or extensions that
  18. they make, so that these may be included in future releases; and (b)
  19. to inform MIT of noteworthy uses of this software.
  20.  
  21. 3. All materials developed as a consequence of the use of this
  22. software shall duly acknowledge such use, in accordance with the usual
  23. standards of acknowledging credit in academic research.
  24.  
  25. 4. MIT has made no warrantee or representation that the operation of
  26. this software will be error-free, and MIT is under no obligation to
  27. provide any services, by way of maintenance, update, or otherwise.
  28.  
  29. 5. In conjunction with products arising from the use of this material,
  30. there shall be no use of the name of the Massachusetts Institute of
  31. Technology nor of any adaptation thereof in any advertising,
  32. promotional, or sales literature without prior written consent from
  33. MIT in each case. */
  34.  
  35. /* Time-Frequency Transforms (pas) */
  36.  
  37. #include "scheme.h"
  38. #include "prims.h"
  39. #include "zones.h"
  40. #include <math.h>
  41. #include "array.h"
  42. #include "image.h"
  43.  
  44. /* SUMMARY
  45.    - pas_cft  (complex data, DIF, split-radix)
  46.    - pas_rft  (real data,    DIT, split-radix) output is conjugate-symmetric
  47.    - pas_csft (cs data,      DIF, split-radix) output is real
  48.    - pas_cft
  49.    - pas_rft2d
  50.    - pas_csft2d
  51.  
  52.  
  53.    Stuff before 4-15-1989
  54.    - C_Array_FFT  (complex data, radix=2, NOT-in-place)
  55.    - CZT (chirp-z-transform) uses the old cft (hence slow).
  56.    - 2d DFT
  57.    */
  58.  
  59. /* The DFT is as defined in Siebert 6003 book,
  60.    i.e.
  61.    forward DFT   =  Negative exponent and division by N
  62.    backward DFT  =  Positive exponent
  63.    (note Seibert forward DFT is Oppenheim backward DFT)
  64.    */
  65.  
  66. /* mathematical constants */
  67. #ifdef PI
  68. #undef PI
  69. #endif
  70. #define PI    3.141592653589793238462643
  71. #define TWOPI 6.283185307179586476925287
  72. #define SQRT_2          1.4142135623730950488
  73. #define ONE_OVER_SQRT_2  .7071067811865475244
  74. /* Abramowitz and Stegun */
  75.  
  76. DEFINE_PRIMITIVE ("PAS-CFT!", Prim_pas_cft, 5, 5, 0)
  77. { long i, length, power, flag;
  78.   REAL *f1,*f2,  *wcos,*w3cos,*w3sin;
  79.   void pas_cft();
  80.   PRIMITIVE_HEADER (5);
  81.   CHECK_ARG (2, ARRAY_P);    /* real part */
  82.   CHECK_ARG (3, ARRAY_P);    /* imag part */
  83.   CHECK_ARG (4, ARRAY_P);    /* twiddle tables, total length = 3*(length/4)  */
  84.   CHECK_ARG (5, FIXNUM_P);    /* (1)=tables precomputed, else recompute */
  85.  
  86.   flag = (arg_integer (1));
  87.   length = ARRAY_LENGTH(ARG_REF(2));
  88.   if (length != (ARRAY_LENGTH(ARG_REF(3)))) error_bad_range_arg(2);
  89.  
  90.   for (power=0, i=length; i>1; power++)
  91.   { if ( (i % 2) == 1) error_bad_range_arg(2);
  92.     i=i/2; }
  93.  
  94.   f1 = ARRAY_CONTENTS(ARG_REF(2));
  95.   f2 = ARRAY_CONTENTS(ARG_REF(3));
  96.   if (f1==f2) error_wrong_type_arg(2);
  97.  
  98.   wcos = ARRAY_CONTENTS(ARG_REF(4)); /* twiddle tables */
  99.   if (ARRAY_LENGTH(ARG_REF(4)) != (3*length/4)) error_bad_range_arg(4);
  100.   w3cos = wcos  + length/4;
  101.   w3sin = w3cos + length/4;
  102.   if ((arg_nonnegative_integer(5)) == 1)
  103.     pas_cft(1, flag, f1,f2, length, power, wcos,w3cos,w3sin);
  104.   else
  105.     pas_cft(0, flag, f1,f2, length, power, wcos,w3cos,w3sin);
  106.   /*        1 means tables are already made
  107.         0 means compute new tables */
  108.  
  109.   PRIMITIVE_RETURN (UNSPECIFIC);
  110. }
  111.  
  112. DEFINE_PRIMITIVE ("PAS-CFT-MAKE-TWIDDLE-TABLES!",
  113.           Prim_pas_cft_make_twiddle_tables, 2, 2, 0)
  114. { long length, power, i;
  115.   REAL  *wcos,*w3cos,*w3sin;
  116.   void pas_cft_make_twiddle_tables_once();
  117.   PRIMITIVE_HEADER (2);
  118.  
  119.   length = arg_nonnegative_integer(1); /* length of cft that we intend to compute */
  120.   CHECK_ARG (2, ARRAY_P);    /*        storage for twiddle tables    */
  121.   if (ARRAY_LENGTH(ARG_REF(2)) != (3*length/4)) error_bad_range_arg(2);
  122.  
  123.   power=0;
  124.   for (power=0, i=length; i>1; power++)
  125.   { if ( (i % 2) == 1) error_bad_range_arg(1);
  126.     i=i/2; }
  127.  
  128.   wcos = ARRAY_CONTENTS(ARG_REF(2)); /* twiddle tables */
  129.   w3cos = wcos  + length/4;
  130.   w3sin = w3cos + length/4;
  131.   pas_cft_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin);
  132.  
  133.   PRIMITIVE_RETURN (UNSPECIFIC);
  134. }
  135.  
  136. /*
  137.   C COMPLEX FOURIER TRANSFORM  (Split-Radix, Decimation-in-frequency)
  138.   C (adapted and optimized from Sorensen,et.al. ASSP-34 no.1 page 152,  February 1986)
  139.   */
  140.  
  141. /* Twiddle Tables for PAS_CFT;
  142.    (tables for forward transform only)
  143.    Inverse transform === forward CFT (without 1/N scaling) followed by time-reversal.
  144.    /
  145.    The tables contain  (2pi/N)*i  for  i=0,1,2,..,N/4
  146.    (except i=0 is ignored, never used)
  147.    /
  148.    Table for wsin[i] is not needed because wsin[i]=wcos[n4-i].
  149.    Table for w3sin[i] is needed however.  The previous relationship does not work for w3sin.
  150.    */
  151.  
  152. /* There are two routines for making twiddle tables:
  153.    a fast one, and a slower one but more precise.
  154.    The differences in speed and accuracy are actually rather small, but anyway.
  155.    Use the slow one for making permanent tables.
  156.    */
  157.  
  158. void
  159. pas_cft_make_twiddle_tables (n,m, wcos,w3cos,w3sin) /* efficient version */
  160.      REAL *wcos, *w3cos, *w3sin;
  161.      long n,m;
  162. { long i, n4;
  163.   double tm;
  164.   REAL costm,sintm;
  165.   n4 = n/4;
  166.   for (i=1; i<n4; i++)        /* start from table entry 1 */
  167.   { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
  168.     wcos[i] = (REAL) cos(tm);
  169.   }
  170.   for (i=1; i<n4; i++)
  171.   { costm = wcos[i];
  172.     sintm = wcos[n4-i];
  173.     w3cos[i] = costm * (1 - 4*sintm*sintm); /* see my notes */
  174.     w3sin[i] = sintm * (4*costm*costm - 1);
  175.   }
  176. }
  177.  
  178. void
  179. pas_cft_make_twiddle_tables_once (n,m, wcos,w3cos,w3sin) /* slow version, more accurate */
  180.      REAL *wcos, *w3cos, *w3sin;
  181.      long n,m;
  182. { long i, n4;
  183.   double tm;
  184.   REAL costm,sintm;
  185.   n4 = n/4;
  186.   for (i=1; i<n4; i++)        /* start from table entry 1 */
  187.   { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
  188.     wcos[i] = (REAL) cos(tm);
  189.     tm = tm * 3.0;        /* this is more precise (in the 16th decimal) than */
  190.     w3cos[i] = (REAL) cos(tm);    /* the more efficient version. (I tested by for/backward) */
  191.     w3sin[i] = (REAL) sin(tm);
  192.   }
  193. }
  194.  
  195. void
  196. pas_cft (tables_ok,flag, x,y,n,m, wcos,w3cos,w3sin)
  197.      REAL *x,*y, *wcos,*w3cos,*w3sin;
  198.      long n,m, flag, tables_ok;
  199. { REAL scale;
  200.   long i,j;
  201.   void pas_cft_make_twiddle_tables();
  202.   void C_Array_Time_Reverse();
  203.   void pas_cft_forward_loop();
  204.  
  205.   if (tables_ok != 1)         /* 1 means = tables already made */
  206.     pas_cft_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
  207.  
  208.   if (flag == 1)        /* forward cft */
  209.   { pas_cft_forward_loop(x,y,n,m, wcos,w3cos,w3sin);
  210.     scale = (REAL) (1.0 / ((double) n));
  211.     for (i=0; i<n; i++)
  212.     { x[i] = x[i]*scale;
  213.       y[i] = y[i]*scale; }}
  214.   else                /* backward cft */
  215.   { for (j=0; j<n; j++) y[j] = (-y[j]); /* conjugate before */
  216.     pas_cft_forward_loop(x,y,n,m, wcos,w3cos,w3sin);
  217.     for (j=0; j<n; j++) y[j] = (-y[j]);    /* conjugate after */
  218.   }
  219. }
  220.  
  221. void
  222. pas_cft_forward_loop (x,y,n,m, wcos,w3cos,w3sin)    /* n >= 4 */
  223.      REAL *x,*y, *wcos,*w3cos,*w3sin;
  224.      long n,m;
  225. { /* REAL  a,a3,e;  no need anymore, use tables */
  226.   REAL    r1,r2,s1,s2,s3,  xt,    cc1,cc3,ss1,ss3;
  227.   long  n1,n2,n4,   i,j,k,    is,id, i0,i1,i2,i3;
  228.   long windex0, windex, windex_n4; /* indices for twiddle tables */
  229.   /********** fortran indices start from 1,... **/
  230.   x = x-1;            /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
  231.   y = y-1;
  232.   /********** fortran indices start from 1,... **/
  233.   /* c */
  234.   /* c-----first M-1 stages of transform */
  235.   /* c */
  236.   windex_n4 = n/4;        /* need for indexing sin via wcos twiddle table */
  237.   n2 = 2*n;
  238.   for (k=1; k<m; k++)        /*  DO 10 K = 1, M-1 */
  239.   { n2 = n2>>1;            /* n2 = n2/2; */
  240.     n4 = n2>>2;            /* n4 = n2/4; */
  241.     /* e = 6.283185307179586476925287 / ((REAL) n2);    no need anymore, use tables */
  242.     /* a = 0.0; */
  243.   {
  244.     /* j=1;  */
  245.     /* The first iteration in the loop "DO 20 J = 1, N4"
  246.        is done specially to save operations involving sin=0, cos=1  */
  247.     /* a = j*e;                         no need anymore, use tables */
  248.     is = 1;            /* is = j; */
  249.     id = 2*n2;
  250.     label40first:
  251.     for (i0=is; i0<n; i0=i0+id) /*  40         DO 30 I0 = IS,N-1,ID */
  252.     { i1 = i0 + n4;
  253.       i2 = i1 + n4;
  254.       i3 = i2 + n4;
  255.       /*    c     */
  256.       r1    = x[i0] - x[i2];
  257.       x[i0] = x[i0] + x[i2];
  258.       r2    = x[i1] - x[i3];
  259.       x[i1] = x[i1] + x[i3];
  260.       s1    = y[i0] - y[i2];
  261.       y[i0] = y[i0] + y[i2];
  262.       s2    = y[i1] - y[i3];
  263.       y[i1] = y[i1] + y[i3];
  264.       /*    c     */
  265.       s3    = r1 - s2;
  266.       r1    = r1 + s2;
  267.       s2    = s1 - r2;        /* original used to be s2 = r2 - s1; */
  268.       r2    = r2 + s1;
  269.       x[i2] =   r1;
  270.       y[i2] =   s2;        /* used to be y[i2] =  (-s2); */
  271.       x[i3] =   s3;
  272.       y[i3] =   r2;
  273.       /* x[i2] =   r1*cc1 + s2*ss1;   used to be, see below
  274.      y[i2] =   s2*cc1 - r1*ss1;   used to be, see below, inside the DO 20 J=1,N4
  275.      x[i3] =   s3*cc3 + r2*ss3;
  276.      y[i3] =   r2*cc3 - s3*ss3; */
  277.     }                /* 30         CONTINUE */
  278.     is = 2*id - n2 + 1;        /* is = 2*id - n2 + j; */
  279.     id = 4*id;
  280.     if (is < n) goto label40first; /* IF (IS.LT.N) GOTO 40 */
  281.   }
  282.     /*  c  */
  283.     windex0 = 1<<(k-1);
  284.     windex  = windex0;
  285.     for (j=2; j<=n4; j++)    /* DO 20 J = 1, N4 */
  286.     {
  287.       /* windex = (j-1)*(1<<(k-1)); -- done with trick to avoid (j-1) and 1<<(k-1) */
  288.       cc1 = wcos[windex];
  289.       ss1 = wcos[windex_n4 - windex]; /* see my notes */
  290.       cc3 = w3cos[windex];
  291.       ss3 = w3sin[windex];    /* sin-from-cos trick does not work here */
  292.       windex = j*windex0;    /* same trick as "a = j*e" */
  293.       /* a3 = 3*a;
  294.      cc1 = cos(a);
  295.      ss1 = sin(a);
  296.      cc3 = cos(a3);
  297.      ss3 = sin(a3);
  298.      a = j*e;*/
  299.       is = j;
  300.       id = 2*n2;
  301.       label40:
  302.       for (i0=is; i0<n; i0=i0+id) /*  40         DO 30 I0 = IS,N-1,ID */
  303.       { i1 = i0 + n4;
  304.     i2 = i1 + n4;
  305.     i3 = i2 + n4;
  306.     /*    c     */
  307.     r1    = x[i0] - x[i2];
  308.     x[i0] = x[i0] + x[i2];
  309.     r2    = x[i1] - x[i3];
  310.     x[i1] = x[i1] + x[i3];
  311.     s1    = y[i0] - y[i2];
  312.     y[i0] = y[i0] + y[i2];
  313.     s2    = y[i1] - y[i3];
  314.     y[i1] = y[i1] + y[i3];
  315.     /*    c     */
  316.     s3    = r1 - s2;
  317.     r1    = r1 + s2;
  318.     s2    = s1 - r2;    /* original used to be s2 = r2 - s1; */
  319.     r2    = r2 + s1;
  320.     x[i2] =   r1*cc1 + s2*ss1; /* used to be x[i2] =   r1*cc1 - s2*ss1;  */
  321.     y[i2] =   s2*cc1 - r1*ss1; /* used to be y[i2] = (-s2*cc1 - r1*ss1); */
  322.     x[i3] =   s3*cc3 + r2*ss3;
  323.     y[i3] =   r2*cc3 - s3*ss3;
  324.       }                /* 30         CONTINUE */
  325.       is = 2*id - n2 + j;
  326.       id = 4*id;
  327.       if (is < n) goto label40; /* IF (IS.LT.N) GOTO 40 */
  328.     }                /* 20      CONTINUE */
  329.   }                /* 10   CONTINUE */
  330.   /* c
  331.      c-----------last-stage, length-2 butterfly ----------------c
  332.      c  */
  333.   is = 1;
  334.   id = 4;
  335.   label50:
  336.   for (i0=is; i0<=n; i0=i0+id)    /* 50   DO 60 I0 = IS, N, ID  */
  337.   { i1    = i0 + 1;
  338.     r1    = x[i0];
  339.     x[i0] = r1 + x[i1];
  340.     x[i1] = r1 - x[i1];
  341.     r1    = y[i0];
  342.     y[i0] = r1 + y[i1];
  343.     y[i1] = r1 - y[i1];
  344.   }                /* 60   CONTINUE */
  345.   is = 2*id - 1;
  346.   id = 4*id;
  347.   if (is < n) goto label50;    /* IF (IS.LT.N) GOTO 50 */
  348.   /*
  349.     c
  350.     c-----------bit-reverse-counter---------------c
  351.     */
  352.   label100:
  353.   j = 1;
  354.   n1 = n - 1;
  355.   for (i=1; i<=n1; i++)        /* DO 104 I = 1, N1 */
  356.   { if (i >= j) goto label101;    /* if (i .ge. j) goto 101 */
  357.     xt = x[j];
  358.     x[j] = x[i];
  359.     x[i] = xt;
  360.     xt = y[j];
  361.     y[j] = y[i];
  362.     y[i] = xt;
  363.     label101: k = n>>1;        /* k = n/2; */
  364.     label102: if (k>=j) goto label103;
  365.     j = j - k;
  366.     k = k>>1;            /* k = k/2; */
  367.     goto label102;
  368.     label103: j = j + k;
  369.   }                /* 104  CONTINUE */
  370.   /* c-------------------------------------*/
  371.   /* c */
  372. }                /* RETURN  END */
  373.  
  374. DEFINE_PRIMITIVE ("PAS-RFT-CSFT!", Prim_pas_rft_csft, 5, 5, 0)
  375. { long i, length, power, flag, ft_type;
  376.   REAL *f1,  *wcos,*w3cos,*w3sin;
  377.   void pas_rft(), pas_csft();
  378.   PRIMITIVE_HEADER (5);
  379.   CHECK_ARG (2, ARRAY_P);    /* Input data (real or cs) */
  380.   CHECK_ARG (3, ARRAY_P);    /* Twiddle tables, total length = 4*(length/8)  */
  381.   CHECK_ARG (4, FIXNUM_P);    /* (1)=tables precomputed, else recompute */
  382.   CHECK_ARG (5, FIXNUM_P);    /* ft_type = 1 or 3
  383.                    1 means compute rft, 3 means compute csft */
  384.   flag = (arg_integer (1));
  385.   f1   = ARRAY_CONTENTS(ARG_REF(2));
  386.   length = ARRAY_LENGTH(ARG_REF(2));
  387.   for (power=0, i=length; i>1; power++)
  388.   { if ( (i % 2) == 1) error_bad_range_arg(2);
  389.     i=i/2; }
  390.  
  391.   wcos = ARRAY_CONTENTS(ARG_REF(3)); /* twiddle tables */
  392.   if (ARRAY_LENGTH(ARG_REF(3)) != (4*length/8)) error_bad_range_arg(3);
  393.   w3cos = wcos + (length/4);
  394.   w3sin = w3cos + (length/8);
  395.  
  396.   ft_type = (arg_nonnegative_integer(5)); /*         rft or csft */
  397.   if (ft_type == 1) {
  398.     if ((arg_nonnegative_integer(4)) == 1)
  399.       pas_rft     (1, flag, f1, length, power, wcos,w3cos,w3sin);
  400.     else pas_rft  (0, flag, f1, length, power, wcos,w3cos,w3sin);
  401.   }
  402.   else if (ft_type == 3) {
  403.     if ((arg_nonnegative_integer(4)) == 1)
  404.       pas_csft    (1, flag, f1, length, power, wcos,w3cos,w3sin);
  405.     else pas_csft (0, flag, f1, length, power, wcos,w3cos,w3sin);
  406.     /*             1 means tables are already made
  407.            0 means compute new tables */
  408.   }
  409.   else error_bad_range_arg(5);
  410.  
  411.   PRIMITIVE_RETURN (UNSPECIFIC);
  412. }
  413.  
  414. DEFINE_PRIMITIVE ("PAS-REALDATA-MAKE-TWIDDLE-TABLES!",
  415.           Prim_pas_realdata_make_twiddle_tables, 2, 2, 0)
  416. { long length, power, i;
  417.   REAL  *wcos,*w3cos,*w3sin;
  418.   void pas_realdata_make_twiddle_tables_once();
  419.   PRIMITIVE_HEADER (2);
  420.  
  421.   length = arg_nonnegative_integer(1); /* length of rft that we intend to compute */
  422.   CHECK_ARG (2, ARRAY_P);    /*        storage for twiddle tables    */
  423.   if (ARRAY_LENGTH(ARG_REF(2)) != (4*length/8)) error_bad_range_arg(2);
  424.  
  425.   power=0;
  426.   for (power=0, i=length; i>1; power++)
  427.   { if ( (i % 2) == 1) error_bad_range_arg(1);
  428.     i=i/2; }
  429.  
  430.   wcos = ARRAY_CONTENTS(ARG_REF(2)); /* twiddle tables */
  431.   w3cos = wcos +  length/4;
  432.   w3sin = w3cos + length/8;
  433.   pas_realdata_make_twiddle_tables_once(length, power, wcos,w3cos,w3sin);
  434.  
  435.   PRIMITIVE_RETURN (UNSPECIFIC);
  436. }
  437.  
  438. /*
  439.   C REAL FOURIER TRANSFORM  (Split-Radix, Decimation-in-time)
  440.   C (adapted from Sorensen,et.al. ASSP-35 no.6 page 849,  October 1986)
  441.   C
  442.   C the output is [Re(0),Re(1),...,Re(n/2), Im(n/2-1),...,Im(1)]
  443.   */
  444.  
  445. /* Twiddle Tables for PAS_RFT and PAS_CSFT
  446.    are identical. -> pas_realdata_make_twiddle_tables
  447.    (but they are indexed differently in each case)
  448.    /
  449.    The tables contain  (2pi/N)*i  where i=0,1,2,..,N/4     for wcos
  450.                   and  (2pi/N)*i  where i=0,1,2,..,N/8     for w3cos w3sin
  451.    (except i=0 is ignored, never used)
  452.    /
  453.    Table for wsin[i] is not needed because wsin[i]=wcos[n4-i].
  454.    Table for w3sin[i] is needed however.  The previous relationship does not work for w3sin.
  455.    /
  456.    Instead of getting sin() from   a wsin[i] table with i=1,..,N/8
  457.    we get it from wcos[n4-i].
  458.    This way we can use a CFT table which goes up to N/4
  459.    for RFT CSFT also. We do so in image-processing (rft2d-csft2d).
  460.    */
  461.  
  462. /* There are two routines for making twiddle tables:
  463.    a fast one, and a slower one but more precise.
  464.    The differences in speed and accuracy are actually rather small, but anyway.
  465.    Use the slow one for making tables that stay around.
  466.    */
  467.  
  468. void pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin)  /* efficient version */
  469.      REAL *wcos, *w3cos, *w3sin;
  470.      long n,m;
  471. { long i, n4, n8;
  472.   double tm;
  473.   REAL costm,sintm;
  474.   n4 = n/4;
  475.   n8 = n/8;
  476.   for (i=1; i<n4; i++)        /* start from table entry 1 */
  477.   { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
  478.     wcos[i] = (REAL) cos(tm);
  479.   }
  480.   for (i=1; i<n8; i++)
  481.   { costm = wcos[i];
  482.     sintm = wcos[n4-i];
  483.     w3cos[i] = costm * (1 - 4*sintm*sintm); /* see my notes */
  484.     w3sin[i] = sintm * (4*costm*costm - 1);
  485.   }
  486. }
  487.  
  488. void pas_realdata_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin) /* slow version, more accurate */
  489.      REAL *wcos, *w3cos, *w3sin;
  490.      long n,m;
  491. { long i, n4, n8;
  492.   double tm;
  493.   REAL costm,sintm;
  494.   n4 = n/4;
  495.   n8 = n/8;
  496.   for (i=1; i<n8; i++)        /* start from table entry 1 */
  497.   { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
  498.     wcos[i] = (REAL) cos(tm);
  499.     tm = tm * 3.0;        /* this is more precise (in the 16th decimal) than */
  500.     w3cos[i] = (REAL) cos(tm);    /* the more efficient version. (I tested by for/backward) */
  501.     w3sin[i] = (REAL) sin(tm);
  502.   }
  503.   for (i=n8; i<n4; i++)
  504.   { tm = 6.283185307179586476925287 * (((double) i) / ((double) n));
  505.     wcos[i] = (REAL) cos(tm);
  506.   }
  507. }
  508.  
  509. void pas_rft(tables_ok,flag, x,n,m, wcos,w3cos,w3sin)
  510.      REAL *x, *wcos,*w3cos,*w3sin;
  511.      long n,m, flag, tables_ok;
  512. { REAL scale;
  513.   long i;
  514.   void pas_realdata_make_twiddle_tables();
  515.   void pas_rft_forward_loop();
  516.  
  517.   if (tables_ok != 1)        /* 1 means = tables already made */
  518.     pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
  519.  
  520.   pas_rft_forward_loop(x,n,m, wcos,w3cos,w3sin);
  521.  
  522.   if (flag == 1)        /* forward rft */
  523.   { scale = (REAL) (1.0 / ((double) n));
  524.     for (i=0; i<n; i++)         x[i] = x[i] * scale; }
  525.   else                /* backward rft */
  526.     for (i=((n/2)+1); i<n; i++) x[i] = (-x[i]); /* time-reverse cs-array */
  527. }
  528.  
  529. /* rft
  530.    forward transform === forward_loop + 1/N scaling
  531.    inverse transform === forward_loop + time-reversal (without 1/N scaling)
  532.    */
  533.  
  534. /* wcos           must be length n/4
  535.    w3cos, w3sin   must be length n/8
  536.    (greater than n/8 is fine also, e.g. use cft tables)
  537.    */
  538.  
  539. void pas_rft_forward_loop(x,n,m, wcos,w3cos,w3sin)
  540.      REAL *x, *wcos,*w3cos,*w3sin;
  541.      long n,m;
  542. { /* REAL   a,a3,e;      no need anymore, use tables */
  543.   REAL   r1, xt,    cc1,cc3,ss1,ss3, t1,t2,t3,t4,t5,t6;
  544.   long n1,n2,n4,n8,  i,j,k,  is,id,   i0,i1,i2,i3,i4,i5,i6,i7,i8;
  545.   long windex0, windex, windex_n4; /* indices for twiddle tables */
  546.   /********** fortran indices start from 1,... **/
  547.   x = x-1;            /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
  548.   /********** fortran indices start from 1,... **/
  549.   /* c */
  550.   windex_n4 = n/4;        /* need for indexing sin via wcos twiddle table */
  551.   /* c
  552.      c-----------bit-reverse-counter---------------c
  553.      */
  554.   label100:
  555.   j = 1;
  556.   n1 = n - 1;
  557.   for (i=1; i<=n1; i++)        /* DO 104 I = 1, N1 */
  558.   { if (i >= j) goto label101;    /* if (i .ge. j) goto 101 */
  559.     xt = x[j];
  560.     x[j] = x[i];
  561.     x[i] = xt;
  562.     label101: k = n>>1;        /* k = n/2; */
  563.     label102: if (k>=j) goto label103;
  564.     j = j - k;
  565.     k = k>>1;            /* k = k/2; */
  566.     goto label102;
  567.     label103: j = j + k;
  568.   }                /* 104  CONTINUE */
  569.   /* c-------------------------------------*/
  570.   /* c */
  571.   /* c  ----length-two-butterflies----------- */
  572.   is = 1;
  573.   id = 4;
  574.   label70:
  575.   for (i0=is; i0<=n; i0=i0+id)  /*  70   DO 60 I0 = IS,N,ID */
  576.   { i1    = i0 + 1;
  577.     r1    = x[i0];
  578.     x[i0] = r1 + x[i1];
  579.     x[i1] = r1 - x[i1];
  580.   }                /* 60   CONTINUE */
  581.   is = 2*id - 1;
  582.   id = 4*id;
  583.   if (is < n) goto label70;    /* IF (IS.LT.N) GOTO 70 */
  584.   /* C
  585.      C -------L-shaped-butterflies-------- */
  586.   n2 = 2;
  587.   for (k=2; k<=m; k++)        /* DO 10 K = 2,M */
  588.   { n2 = n2 * 2;
  589.     n4 = n2>>2;            /* n4 = n2/4; */
  590.     n8 = n2>>3;            /* n8 = n2/8; */
  591.     /* e = 6.283185307179586476925287 / ((REAL) n2);   no need anymore, use tables */
  592.     is = 0;
  593.     id = n2 * 2;
  594.     label40:
  595.     for (i=is; i<n; i=i+id)    /* 40      DO 38 I = IS,N-1,ID */
  596.     { i1 = i + 1;
  597.       i2 = i1 + n4;
  598.       i3 = i2 + n4;
  599.       i4 = i3 + n4;
  600.       t1 = x[i4] + x[i3];
  601.       x[i4] = x[i4] - x[i3];
  602.       x[i3] = x[i1] - t1;
  603.       x[i1] = x[i1] + t1;
  604.       if (n4 == 1) goto label38; /* IF (N4.EQ.1) GOTO 38 */
  605.       i1 = i1 + n8;
  606.       i2 = i2 + n8;
  607.       i3 = i3 + n8;
  608.       i4 = i4 + n8;
  609.       /* t1    = (x[i3] + x[i4]) / sqrt(2.0); -- this is more precise, it uses extended
  610.      t2    = (x[i3] - x[i4]) / sqrt(2.0); -- precision inside 68881, but slower */
  611.       t1    = (x[i3] + x[i4]) * ONE_OVER_SQRT_2;
  612.       t2    = (x[i3] - x[i4]) * ONE_OVER_SQRT_2;
  613.       x[i4] = x[i2] - t1;
  614.       x[i3] = -x[i2] - t1;
  615.       x[i2] = x[i1] - t2;
  616.       x[i1] = x[i1] + t2;
  617.       label38:            /* 38      CONTINUE */
  618.       ;
  619.     }
  620.     is = 2*id - n2;
  621.     id = 4*id;
  622.     if (is < n) goto label40;    /* IF (IS.LT.N) GOTO 40 */
  623.     /* a = e; */
  624.     windex0 = 1<<(m-k);
  625.     windex  = windex0;
  626.     for (j=2; j<=n8; j++)    /* DO 32 J = 2,N8 */
  627.     {
  628.       /* windex = (j-1)*(1<<(m-k)); -- done with trick to avoid (j-1) and 1<<(m-k) */
  629.       cc1 = wcos[windex];
  630.       ss1 = wcos[windex_n4 - windex]; /* sin-from-cos trick: see my notes */
  631.       cc3 = w3cos[windex];
  632.       ss3 = w3sin[windex];    /* sin-from-cos trick does not work here */
  633.       windex = j*windex0;    /* same trick as "a = j*e" */
  634.       /* a3 = 3*a;
  635.      cc1 = cos(a);
  636.      ss1 = sin(a);
  637.      cc3 = cos(a3);
  638.      ss3 = sin(a3);
  639.      a = j*e;*/
  640.       is = 0;
  641.       id = 2*n2;
  642.       label36:            /*  36         DO 30 I = IS,N-1,ID */
  643.       for (i=is; i<n; i=i+id)
  644.       { i1 = i + j;
  645.     i2 = i1 + n4;
  646.     i3 = i2 + n4;
  647.     i4 = i3 + n4;
  648.     i5 = i  + n4 - j + 2;
  649.     i6 = i5 + n4;
  650.     i7 = i6 + n4;
  651.     i8 = i7 + n4;
  652.     t1 = x[i3]*cc1 + x[i7]*ss1;
  653.     t2 = x[i7]*cc1 - x[i3]*ss1;
  654.     t3 = x[i4]*cc3 + x[i8]*ss3;
  655.     t4 = x[i8]*cc3 - x[i4]*ss3;
  656.     t5 = t1 + t3;
  657.     t6 = t2 + t4;
  658.     t3 = t3 - t1;        /* t3 = t1 - t3; */
  659.     t4 = t2 - t4;
  660.     x[i8] = x[i6] + t6;
  661.     x[i3] = t6    - x[i6];
  662.     x[i4] = x[i2] + t3;    /* x[i4] = x[i2] - t3; */
  663.     x[i7] = t3 - x[i2];    /* x[i7] = -x[i2] - t3; */
  664.     x[i6] = x[i1] - t5;
  665.     x[i1] = x[i1] + t5;
  666.     x[i2] = x[i5] + t4;
  667.     x[i5] = x[i5] - t4;
  668.       }                /* 30         CONTINUE */
  669.       is = 2*id - n2;
  670.       id = 4*id;
  671.       if (is < n) goto label36; /* IF (IS.LT.N) GOTO 36 */
  672.     }                /* 32      CONTINUE */
  673.   }                /* 10   CONTINUE */
  674. }                /* RETURN  END */
  675.  
  676.  
  677. /*
  678.   C CONJUGATE SYMMETRIC FOURIER TRANSFORM  (Split-Radix, Decimation-in-time)
  679.   C (adapted from Sorensen,et.al. ASSP-35 no.6 page 849,  October 1986)
  680.   C
  681.   C input is [Re(0),Re(1),...,Re(n/2), Im(n/2-1),...,Im(1)]
  682.   C output is real
  683.   */
  684.  
  685. /* twiddle tables identical with rft
  686.    for comments see rft */
  687.  
  688. void pas_csft(tables_ok,flag, x,n,m, wcos,w3cos,w3sin)
  689.      REAL *x, *wcos,*w3cos,*w3sin;
  690.      long n,m, flag, tables_ok;
  691. { REAL scale;
  692.   long i,n2;
  693.   void pas_realdata_make_twiddle_tables();
  694.   void C_Array_Time_Reverse();
  695.   void pas_csft_backward_loop();
  696.  
  697.   if (tables_ok != 1)        /* 1 means = tables already made */
  698.     pas_realdata_make_twiddle_tables(n,m, wcos,w3cos,w3sin);
  699.  
  700.   if (flag == 1)        /* forward csft */
  701.   { n2 = n/2;
  702.     scale = (REAL) (1.0 / ((double) n));
  703.     for (i=0; i<=n2; i++)   x[i] = x[i]*scale;
  704.     scale = (-scale);
  705.     for (i=n2+1; i<n; i++)  x[i] = x[i]*scale; /* scale and conjugate cs-array */
  706.     pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin);
  707.   }
  708.   else                /* backward csft */
  709.     pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin);
  710. }
  711.  
  712. /* csft
  713.    forward transform === backward_loop + 1/N scaling + time-reversal
  714.    inverse transform === backward_loop
  715.    */
  716.  
  717. /* wcos           must be length n/4
  718.    w3cos, w3sin   must be length n/8
  719.    (greater than n/8 is fine also, e.g. use cft tables)
  720.    */
  721.  
  722. void pas_csft_backward_loop(x,n,m, wcos,w3cos,w3sin)
  723.      REAL *x, *wcos,*w3cos,*w3sin;
  724.      long n,m;
  725. { /* REAL   a,a3,e;     no need anymore, use tables */
  726.   REAL   r1, xt,    cc1,cc3,ss1,ss3, t1,t2,t3,t4,t5;
  727.   long n1,n2,n4,n8,  i,j,k,  is,id,   i0,i1,i2,i3,i4,i5,i6,i7,i8;
  728.   long windex0, windex, windex_n4; /* indices for twiddle tables */
  729.   /********** fortran indices start from 1,... **/
  730.   x = x-1;            /* TRICK---- x(0) is now illegal, but x(1) and x(n) are valid */
  731.   /********** fortran indices start from 1,... **/
  732.   /* c */
  733.   windex_n4 = n/4;        /* need for indexing sin via wcos twiddle table */
  734.   /* c */
  735.   /* c */
  736.   /* c
  737.      c -------L-shaped-butterflies-------- */
  738.   n2 = 2*n;
  739.   for (k=1; k<m; k++)        /* do 10 k = 1,m-1 */
  740.   { is = 0;
  741.     id = n2;
  742.     n2 = n2>>1;            /* n2 = n2/2; */
  743.     n4 = n2>>2;            /* n4 = n2/4; */
  744.     n8 = n4>>1;            /* n8 = n4/2; */
  745.     /* e = 6.283185307179586476925287 / ((REAL) n2);  no need anymore, use tables */
  746.     label17:
  747.     for (i=is; i<n; i=i+id)    /*  17      do 15 i = is,(n-1),id */
  748.     { i1 = i + 1;
  749.       i2 = i1 + n4;
  750.       i3 = i2 + n4;
  751.       i4 = i3 + n4;
  752.       t1    = x[i1] - x[i3];
  753.       x[i1] = x[i1] + x[i3];
  754.       x[i2] = 2*x[i2];
  755.       x[i3] = t1  - 2*x[i4];
  756.       x[i4] = t1  + 2*x[i4];
  757.       if (n4 == 1) goto label15; /* if (n4.eq.1) goto 15 */
  758.       i1 = i1 + n8;
  759.       i2 = i2 + n8;
  760.       i3 = i3 + n8;
  761.       i4 = i4 + n8;
  762.       t1    =   (x[i2] - x[i1])  * ONE_OVER_SQRT_2;
  763.       t2    = (-(x[i4] + x[i3])) * ONE_OVER_SQRT_2;
  764.       /* t1    = (x[i2] - x[i1])/sqrt(2.0);
  765.      t2    = (x[i4] + x[i3])/sqrt(2.0); */
  766.       x[i1] = x[i1] + x[i2];
  767.       x[i2] = x[i4] - x[i3];
  768.       x[i3] = 2 * (t2-t1);    /* x[i3] = 2 * (-t2-t1); */
  769.       x[i4] = 2 * (t2+t1);    /* x[i4] = 2 * (-t2+t1); */
  770.       label15:
  771.       ;
  772.     }                /* 15      continue */
  773.     is = 2*id - n2;
  774.     id = 4*id;
  775.     if (is < (n-1)) goto label17; /* if (is.lt.(n-1)) goto 17 */
  776.     /* a = e; */
  777.     windex0 = 1<<(k-1);        /* see my notes */
  778.     windex  = windex0;
  779.     for (j=2; j<=n8; j++)    /* do 20 j=2,n8 */
  780.     {
  781.       /* windex = (j-1)*(1<<(k-1)); -- done with trick to avoid (j-1) and 1<<(k-1) */
  782.       cc1 = wcos[windex];
  783.       ss1 = wcos[windex_n4 - windex]; /* sin-from-cos trick: see my notes */
  784.       cc3 = w3cos[windex];
  785.       ss3 = w3sin[windex];    /* sin-from-cos trick does not work here */
  786.       windex = j*windex0;    /* same trick as "a = j*e" */
  787.       /* a3 = 3*a;
  788.      cc1 = cos(a);
  789.      ss1 = sin(a);
  790.      cc3 = cos(a3);
  791.      ss3 = sin(a3);
  792.      a = j*e; */
  793.       is = 0;
  794.       id = 2*n2;
  795.       label40:
  796.       for (i=is; i<n; i=i+id)    /* 40         do 30 i = is,(n-1),id */
  797.       { i1 = i + j;
  798.     i2 = i1 + n4;
  799.     i3 = i2 + n4;
  800.     i4 = i3 + n4;
  801.     i5 = i  + n4 - j + 2;
  802.     i6 = i5 + n4;
  803.     i7 = i6 + n4;
  804.     i8 = i7 + n4;
  805.     t1    = x[i1] - x[i6];
  806.     x[i1] = x[i1] + x[i6];
  807.     t2    = x[i5] - x[i2];
  808.     x[i5] = x[i5] + x[i2];
  809.     t3    = x[i8] + x[i3];
  810.     x[i6] = x[i8] - x[i3];
  811.     t4    = x[i4] + x[i7];
  812.     x[i2] = x[i4] - x[i7];
  813.     t5 = t1 - t4;
  814.     t1 = t1 + t4;
  815.     t4 = t2 - t3;
  816.     t2 = t2 + t3;
  817.     x[i3] = t5*cc1    + t4*ss1;
  818.     x[i7] = t5*ss1    - t4*cc1;     /* x[i7] = (-t4*cc1) + t5*ss1; */
  819.     x[i4] = t1*cc3    - t2*ss3;
  820.     x[i8] = t2*cc3    + t1*ss3;
  821.       }                /* 30         continue */
  822.       is = 2*id - n2;
  823.       id = 4*id;
  824.       if (is < (n-1)) goto label40; /* if (is.lt.(n-1)) goto 40 */
  825.     }                /*  20      continue */
  826.   }                /* 10   continue */
  827.   /* c */
  828.   /* c  ----length-two-butterflies----------- */
  829.   is = 1;
  830.   id = 4;
  831.   label70:
  832.   for (i0=is; i0<=n; i0=i0+id)  /*  70   DO 60 I0 = IS,N,ID */
  833.   { i1    = i0 + 1;
  834.     r1    = x[i0];
  835.     x[i0] = r1 + x[i1];
  836.     x[i1] = r1 - x[i1];
  837.   }                /* 60   CONTINUE */
  838.   is = 2*id - 1;
  839.   id = 4*id;
  840.   if (is < n) goto label70;    /* IF (IS.LT.N) GOTO 70 */
  841.   /* c */
  842.   /* c-----------bit-reverse-counter---------------c */
  843.   label100:
  844.   j = 1;
  845.   n1 = n - 1;
  846.   for (i=1; i<=n1; i++)        /* DO 104 I = 1, N1 */
  847.   { if (i >= j) goto label101;    /* if (i .ge. j) goto 101 */
  848.     xt = x[j];
  849.     x[j] = x[i];
  850.     x[i] = xt;
  851.     label101: k = n>>1;        /* k = n/2; */
  852.     label102: if (k>=j) goto label103;
  853.     j = j - k;
  854.     k = k>>1;            /* k = k/2; */
  855.     goto label102;
  856.     label103: j = j + k;
  857.   }                /* 104  CONTINUE */
  858.   /* c */
  859. }                /* RETURN  END */
  860.  
  861.  
  862.  
  863.  
  864. /* Image processing     only for square images
  865.    (old stuff handles non-square but is slow)
  866.    For 2d FTs  precomputed tables or not, make almost no difference in total time.
  867.    */
  868.  
  869. DEFINE_PRIMITIVE ("PAS-CFT2D!", Prim_pas_cft2d, 5,5, 0)
  870. { long i, length, power, flag, rows,rowpower;
  871.   REAL *f1,*f2,  *wcos,*w3cos,*w3sin;
  872.   void pas_cft2d();
  873.   PRIMITIVE_HEADER (5);
  874.   CHECK_ARG (2, ARRAY_P);    /* real part */
  875.   CHECK_ARG (3, ARRAY_P);    /* imag part */
  876.   CHECK_ARG (4, ARRAY_P);    /* twiddle tables, length = 3*(rows/4)  */
  877.  
  878.   flag = (arg_integer (1));
  879.   length = ARRAY_LENGTH(ARG_REF(2));
  880.   if (length != (ARRAY_LENGTH(ARG_REF(3)))) error_bad_range_arg(2);
  881.  
  882.   for (power=0, i=length; i>1; power++)    /*         length must be power of 2 */
  883.   { if ( (i % 2) == 1) error_bad_range_arg(2);
  884.     i=i/2; }
  885.  
  886.   if ((power % 2) == 1) error_bad_range_arg(2);
  887.   rowpower = (power/2);
  888.   rows = (1<<rowpower);        /*                 square image */
  889.  
  890.   f1 = ARRAY_CONTENTS(ARG_REF(2));
  891.   f2 = ARRAY_CONTENTS(ARG_REF(3));
  892.   if (f1==f2) error_wrong_type_arg(2);
  893.  
  894.   wcos = ARRAY_CONTENTS(ARG_REF(4)); /* twiddle tables */
  895.   if (ARRAY_LENGTH(ARG_REF(4)) != (3*rows/4)) error_bad_range_arg(4);
  896.   w3cos = wcos   +   rows/4;
  897.   w3sin = w3cos  +   rows/4;
  898.   if ((arg_nonnegative_integer(5)) == 1)
  899.     pas_cft2d(1, flag, f1,f2, rows, rowpower, wcos,w3cos,w3sin);
  900.   else
  901.     pas_cft2d(0, flag, f1,f2, rows, rowpower, wcos,w3cos,w3sin);
  902.   /*          1 means tables are already made
  903.           0 means compute new tables */
  904.  
  905.   PRIMITIVE_RETURN (UNSPECIFIC);
  906. }
  907.  
  908. /* pas_cft2d
  909.    n =                rows of square image, rows is power of 2
  910.    m =                rowpower
  911.    Scaling (1/n) is done all-at-once at the end.
  912.    Time-Reversing is done intermediately, it is more efficient.
  913.    */
  914. void pas_cft2d(tables_ok,flag, x,y,n,m, wcos,w3cos,w3sin)
  915.      REAL *x,*y, *wcos,*w3cos,*w3sin;
  916.      long n,m, flag, tables_ok;
  917. { REAL scale, *xrow,*yrow;
  918.   long i,j, rows,cols, total_length;
  919.   void pas_cft_make_twiddle_tables_once();
  920.   void C_Array_Time_Reverse();
  921.   void pas_cft_forward_loop();
  922.  
  923.   if (tables_ok != 1)         /* 1 means = tables already made */
  924.     pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
  925.  
  926.   rows = n;
  927.   cols = rows;            /* square image */
  928.   total_length = rows*rows;
  929.  
  930.   if (flag != 1)        /* backward transform */
  931.     for (i=0; i<total_length; i++) y[i] = (-y[i]); /* conjugate before */
  932.  
  933.   xrow = x; yrow = y;        /* ROW-WISE */
  934.   for (i=0; i<rows; i++)    /* forward or backward */
  935.   { pas_cft_forward_loop( xrow,yrow, n,m, wcos,w3cos,w3sin);
  936.     xrow = xrow + cols;
  937.     yrow = yrow + cols; }
  938.  
  939.   Image_Fast_Transpose(x, rows); /* COLUMN-WISE */
  940.   Image_Fast_Transpose(y, rows);
  941.   xrow = x; yrow = y;
  942.   for (i=0; i<rows; i++)    /* forward or backward */
  943.   { pas_cft_forward_loop( xrow,yrow, n,m, wcos,w3cos,w3sin);
  944.     xrow = xrow + cols;
  945.     yrow = yrow + cols; }
  946.  
  947.   Image_Fast_Transpose(x, rows);
  948.   Image_Fast_Transpose(y, rows);
  949.  
  950.   if (flag == 1)        /* forward : scale */
  951.   { scale = (REAL) (1.0 / ((double) total_length));
  952.     for (i=0; i<total_length; i++)
  953.     { x[i] = x[i]*scale;
  954.       y[i] = y[i]*scale; }}
  955.   else                /* backward : conjugate after */
  956.     for (i=0; i<total_length; i++) y[i] = (-y[i]);
  957. }
  958.  
  959.  
  960. DEFINE_PRIMITIVE ("PAS-RFT2D-CSFT2D!", Prim_pas_rft2d_csft2d, 5,5, 0)
  961. { long i, length, power, flag, ft_type, rows,rowpower;
  962.   REAL *f1,  *wcos,*w3cos,*w3sin;
  963.   void pas_rft2d(), pas_csft2d();
  964.   PRIMITIVE_HEADER (5);
  965.   CHECK_ARG (2, ARRAY_P);    /* Input data (real or cs) */
  966.   CHECK_ARG (3, ARRAY_P);    /* CFT twiddle tables, length = 3*(rows/4)  */
  967.   CHECK_ARG (4, FIXNUM_P);    /* (1)=tables precomputed, else recompute */
  968.   flag = (arg_integer (1));
  969.   f1 = ARRAY_CONTENTS(ARG_REF(2));
  970.   length = ARRAY_LENGTH(ARG_REF(2));
  971.   for (power=0, i=length; i>1; power++)    /* length must be power of 2 */
  972.   { if ( (i % 2) == 1) error_bad_range_arg(2);
  973.     i=i/2; }
  974.  
  975.   if ((power % 2) == 1) error_bad_range_arg(2);
  976.   rowpower = (power/2);
  977.   rows = (1<<rowpower);        /*                 square image */
  978.  
  979.   wcos = ARRAY_CONTENTS(ARG_REF(3)); /* CFT twiddle tables */
  980.   if (ARRAY_LENGTH(ARG_REF(3)) != (3*rows/4)) error_bad_range_arg(3);
  981.   w3cos = wcos  + rows/4;
  982.   w3sin = w3cos + rows/4;
  983.  
  984.   ft_type = (arg_nonnegative_integer(5)); /*          rft2d or csft2d */
  985.   if (ft_type == 1) {
  986.     if ((arg_nonnegative_integer(4)) == 1)
  987.       pas_rft2d     (1, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
  988.     else pas_rft2d  (0, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
  989.   }
  990.   else if (ft_type == 3) {
  991.     if ((arg_nonnegative_integer(4)) == 1)
  992.       pas_csft2d    (1, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
  993.     else pas_csft2d (0, flag, f1, rows, rowpower, wcos,w3cos,w3sin);
  994.     /*               1 means tables are already made
  995.              0 means compute new tables */
  996.   }
  997.   else  error_bad_range_arg(5);
  998.  
  999.   PRIMITIVE_RETURN (UNSPECIFIC);
  1000. }
  1001.  
  1002. /* c                             RFT2D      CSFT2D
  1003.    The frequencies are scrabled wrt  what cft2d (and the old image-fft) give.
  1004.    See   cs-image-magnitude  and  cs-image-real    which unscrable automatically.
  1005.    c
  1006.    c Implementation notes:
  1007.    c   conjugate in one domain         is         reverse and conjugate in other
  1008.    c   reverse   in one domain         is         reverse in other
  1009.    c
  1010.    c    reverse cs-array     which is identical to      conjugate cs-array  (same domain)
  1011.    c                                   is         reverse in other domain
  1012.    c
  1013.    c conjugate cs-array before csft    is-better-than    reverse real-array afterwards
  1014.    c
  1015.    c
  1016.    c    rft2d-csft2d  use 1d-cft tables    to compute rft
  1017.    c    cft tables are simply larger than realdata tables.
  1018.    */
  1019.  
  1020. /* pas_rft2d
  1021.    n =                rows of square image, rows is power of 2
  1022.    m =                rowpower
  1023.    Scaling (1/n) is done all-at-once at the end.
  1024.    Time-Reversing is done intermediately, it is more efficient.
  1025.    */
  1026. void pas_rft2d(tables_ok,flag, x, n,m, wcos,w3cos,w3sin)
  1027.      REAL *x, *wcos,*w3cos,*w3sin;
  1028.      long n,m, flag, tables_ok;
  1029. { REAL scale, *xrow,*yrow;
  1030.   long i,j, rows,cols, total_length, n2;
  1031.   void pas_cft_make_twiddle_tables_once();
  1032.   void C_Array_Time_Reverse();
  1033.   void pas_rft_forward_loop(), pas_cft_forward_loop();
  1034.  
  1035.   if (tables_ok != 1)         /* 1 means = tables already made */
  1036.     pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
  1037.  
  1038.   rows = n;
  1039.   cols = rows;            /* square image */
  1040.   n2   = n/2;
  1041.   total_length = rows*rows;
  1042.  
  1043.   xrow = x;            /*                First ROW-WISE */
  1044.   if (flag == 1)        /* forward transform */
  1045.     for (i=0; i<rows; i++)
  1046.     { pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1047.       xrow = xrow + cols; }
  1048.   else                /* backward transform */
  1049.     for (i=0; i<rows; i++)
  1050.     { pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1051.       for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
  1052.       xrow = xrow + cols; }
  1053.  
  1054.   Image_Fast_Transpose(x, rows); /* COLUMN-WISE */
  1055.  
  1056.   /*      TREAT specially rows 0 and n2,
  1057.       they are real and go into cs-arrays */
  1058.   if (flag == 1)        /* forward transform */
  1059.   { xrow =          x + 0      ;
  1060.     pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1061.     xrow =          x + n2*cols;
  1062.     pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin); }
  1063.   else                /* backward transform */
  1064.   { xrow =          x + 0      ;
  1065.     pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1066.     for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
  1067.     xrow =          x + n2*cols;
  1068.     pas_rft_forward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1069.     for (j=n2+1; j<n; j++) xrow[j] = (-xrow[j]); /* time-reverse cs-array */
  1070.   }
  1071.  
  1072.   /*     TREAT the rest of the rows with CFT
  1073.    */
  1074.   if (flag != 1)        /* backward : conjugate before */
  1075.     for (i=(n2+1)*cols; i<total_length; i++)    x[i] = (-x[i]);
  1076.  
  1077.   xrow = x + cols;        /* real part */
  1078.   yrow = x + (rows-1)*cols;    /* imag part */
  1079.   for (i=1; i<n2; i++)        /* forward or backward transform */
  1080.   { pas_cft_forward_loop(xrow,yrow, n,m, wcos,w3cos,w3sin);
  1081.     xrow = xrow + cols;
  1082.     yrow = yrow - cols; }
  1083.   /*    DO NOT TRANSPOSE BACK, leave real-imag in horizontal rows, save.
  1084.    */
  1085.  
  1086.   if (flag == 1)        /* forward : scale */
  1087.   { scale = (REAL) (1.0 / ((double) total_length));
  1088.     for (i=0; i<total_length; i++)
  1089.       x[i] = x[i]*scale; }
  1090.   else                /* backward : conjugate after */
  1091.     for (i=(n2+1)*cols; i<total_length; i++)
  1092.       x[i] = (-x[i]);
  1093. }
  1094.  
  1095.  
  1096. /* pas_csft2d
  1097.    n =                rows of square image, rows is power of 2
  1098.    m =                rowpower
  1099.    Scaling (1/n) is done all-at-once at the end.
  1100.    Time-Reversing is done intermediately, it is more efficient.
  1101.    */
  1102. void pas_csft2d(tables_ok,flag, x, n,m, wcos,w3cos,w3sin)
  1103.      REAL *x, *wcos,*w3cos,*w3sin;
  1104.      long n,m, flag, tables_ok;
  1105. { REAL scale, *xrow,*yrow;
  1106.   long i,j, rows,cols, total_length, n2;
  1107.   void pas_cft_make_twiddle_tables_once();
  1108.   void C_Array_Time_Reverse();
  1109.   void pas_csft_backward_loop(), pas_cft_forward_loop();
  1110.  
  1111.   if (tables_ok != 1)         /* 1 means = tables already made */
  1112.     pas_cft_make_twiddle_tables_once(n,m, wcos,w3cos,w3sin);
  1113.  
  1114.   rows = n;
  1115.   cols = rows;            /* square image */
  1116.   n2   = n/2;
  1117.   total_length = rows*rows;
  1118.  
  1119.   /*                                     First  ROW-WISE */
  1120.  
  1121.   /*      TREAT SPECIALLY ROWS 0 and n2,   they are cs-arrays and they go into real */
  1122.   if (flag == 1)        /* forward transform */
  1123.   { xrow =          x + 0      ;
  1124.     for (j=n2+1; j<n; j++)  xrow[j]=(-xrow[j]); /* conjugate before */
  1125.     pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1126.     xrow =          x + n2*cols;
  1127.     for (j=n2+1; j<n; j++)  xrow[j]=(-xrow[j]); /* conjugate before */
  1128.     pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin); }
  1129.   else                /* backward transform */
  1130.   { xrow =          x + 0      ;
  1131.     pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1132.     xrow =          x + n2*cols;
  1133.     pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin); }
  1134.  
  1135.   /*     TREAT the rest of the rows with CFT
  1136.    */
  1137.   if (flag != 1)        /* backward : conjugate before */
  1138.     for (i=(n2+1)*cols; i<total_length; i++)    x[i] = (-x[i]);
  1139.  
  1140.   xrow = x + cols;        /* real part */
  1141.   yrow = x + (rows-1)*cols;    /* imag part */
  1142.   for (i=1; i<n2; i++)        /* forward or backward transform */
  1143.   { pas_cft_forward_loop(xrow,yrow, n,m, wcos,w3cos,w3sin);
  1144.     xrow = xrow + cols;
  1145.     yrow = yrow - cols; }
  1146.  
  1147.   if (flag != 1)        /* backward : conjugate after */
  1148.     for (i=(n2+1)*cols; i<total_length; i++)    x[i] = (-x[i]);
  1149.  
  1150.   Image_Fast_Transpose(x, rows);
  1151.   /*                                Second   COLUMN-WISE
  1152.                     Everything should be cs-arrays now */
  1153.  
  1154.   xrow = x;
  1155.   if (flag == 1)        /* forward transform */
  1156.     for (i=0; i<rows; i++)
  1157.     { for (j=n2+1; j<n; j++)  xrow[j]=(-xrow[j]); /* conjugate before */
  1158.       pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1159.       xrow = xrow + cols; }
  1160.   else                /* backward transform */
  1161.     for (i=0; i<rows; i++)
  1162.     { pas_csft_backward_loop(xrow, n,m, wcos,w3cos,w3sin);
  1163.       xrow = xrow + cols; }
  1164.  
  1165.   if (flag == 1)        /* forward : scale */
  1166.   { scale = (REAL) (1.0 / ((double) total_length));
  1167.     for (i=0; i<total_length; i++)
  1168.       x[i] = x[i] * scale; }
  1169. }
  1170.  
  1171.  
  1172.  
  1173. /* STUFF BEFORE 4-15-1989
  1174.  */
  1175.  
  1176. void Make_FFT_Tables(w1, w2, n, flag)
  1177.      REAL *w1, *w2; long n, flag;         /* n  = length of data */
  1178. { long m, n2=n/2;                         /* n2 = length of w1,w2 */
  1179.   double tm;
  1180.   if (flag==1)            /* FORWARD FFT */
  1181.     for (m=0; m<n2; m++) {
  1182.       tm = TWOPI * ((double) m) / ((double) n);
  1183.       w1[m] = (REAL) cos(tm);
  1184.       w2[m] = (REAL) - sin(tm); }
  1185.   else
  1186.     for (m=0; m<n2; m++) {
  1187.       tm = TWOPI * ((double) m) / ((double) n);
  1188.       w1[m] = (REAL) cos(tm);
  1189.       w2[m] = (REAL) sin(tm); }
  1190. }
  1191.  
  1192. #define mult(pf1, pf2, pg1, pg2, w1, w2)            \
  1193.     {  long x, y, p2, p3, p4, p5, p6, p7;           \
  1194.        REAL tmp1, tmp2;                             \
  1195.        a = a / 2;                                   \
  1196.        p2 = - a;                                    \
  1197.        p3 = 0;                                      \
  1198.        for ( x = 1; x <= n2; x = x + a ) {          \
  1199.      p2 = p2 + a;                               \
  1200.      for( y = 1; y <= a; ++y ) {                \
  1201.        ++p3;                                    \
  1202.        p4 = p2 + 1;                             \
  1203.        p5 = p2 + p3;                            \
  1204.        p5 = ((p5-1) % n) + 1;                   \
  1205.        p6 = p5 + a;                             \
  1206.        tmp1 =  w1[p4-1] * pf1[p6-1]             \
  1207.              - w2[p4-1] * pf2[p6-1];            \
  1208.        tmp2 =  w1[p4-1] * pf2[p6-1]             \
  1209.                  + w2[p4-1] * pf1[p6-1];            \
  1210.        pg1[p3-1] = pf1[p5-1] + tmp1;            \
  1211.        pg2[p3-1] = pf2[p5-1] + tmp2;            \
  1212.        p7 = p3 + n2;                            \
  1213.        pg1[p7-1] = pf1[p5-1] - tmp1;            \
  1214.        pg2[p7-1] = pf2[p5-1] - tmp2;            \
  1215.      }                                          \
  1216.        }                                            \
  1217. }
  1218.  
  1219. /* n = length of input data f1,f2;
  1220.    power =  log2(n),
  1221.    g1,g2  are intermediate arrays of length n,
  1222.    w1,w2 point to FFT tables (twiddle factors),
  1223.    flag 1 for forward FFT, else inverse.
  1224.    */
  1225. /* The arrays w1,w2 are half the size of f1,f2,g1,g2.
  1226.    f1,f2 contain the real and imaginary parts of the signal.
  1227.    The answer is left in f1, f2.
  1228.    */
  1229.  
  1230. /* C_Array_FFT
  1231.    complex data, radix=2, not-in-place.
  1232.    (adapted from an fft program I got from Yekta)
  1233.    */
  1234. void C_Array_FFT(flag, power, n, f1, f2, g1,g2,w1,w2)
  1235.      long flag, power, n; REAL f1[], f2[], g1[], g2[], w1[], w2[];
  1236. { long n2=n>>1, a;
  1237.   long  i, l, m;
  1238.   REAL tm;
  1239.   a = n;            /* initially equal to length of data */
  1240.  
  1241.   for (m=0; m<n; m++) { g1[m] = f1[m]; g2[m] = f2[m]; }
  1242.   Make_FFT_Tables(w1,w2,n, flag);
  1243.  
  1244.   if ((power % 2) == 1) l = 2; else l = 1; /* even,odd power */
  1245.   for ( i = l; i <= power ; i = i + 2 ) {
  1246.     mult(g1,g2,f1,f2,w1,w2);
  1247.     mult(f1,f2,g1,g2,w1,w2); }
  1248.  
  1249.   if (flag==1) {        /* FORWARD FFT */
  1250.     tm = 1. / ((REAL) n);    /* normalizing factor */
  1251.     if (l==1)            /* even power */
  1252.       for (m=0; m<n; m++) { f1[m] = tm * g1[m];    f2[m] = tm * g2[m]; }
  1253.     else {            /* odd power ==> do one more mult */
  1254.       mult(g1,g2,f1,f2,w1,w2);    /* f1 and f2 contain the result now */
  1255.       for (m=0; m<n; m++) { f1[m] = tm * f1[m]; f2[m] = tm * f2[m]; }}
  1256.   }
  1257.   else {            /* BACKWARD FFT */
  1258.     if (l==1)            /* even power */
  1259.       for (m=0; m<n; m++) { f1[m] = g1[m]; f2[m] = g2[m]; }
  1260.     else            /* odd power ==> do one more mult */
  1261.       mult(g1,g2,f1,f2,w1,w2);    /* f1 and f2 contain the result now */
  1262.   }
  1263. }
  1264.  
  1265. void C_Array_FFT_With_Given_Tables(flag, power, n, f1, f2, g1,g2,w1,w2)
  1266.      long flag, power, n; REAL f1[], f2[], g1[], g2[], w1[], w2[];
  1267. { long n2=n>>1, a;
  1268.   long  i, l, m;
  1269.   REAL tm;
  1270.   a = n;            /* initially equal to length */
  1271.  
  1272.   for (m=0; m<n; m++) { g1[m] = f1[m];  g2[m] = f2[m]; }
  1273.  
  1274.   if ((power % 2) == 1) l = 2; else l = 1; /* even,odd power */
  1275.   for ( i = l; i <= power ; i = i + 2 ) {
  1276.     mult(g1,g2,f1,f2,w1,w2);
  1277.     mult(f1,f2,g1,g2,w1,w2); }
  1278.  
  1279.   if (flag==1) {        /* FORWARD FFT */
  1280.     tm = 1. / ((REAL) n);    /* normalizing factor */
  1281.     if (l==1)            /* even power */
  1282.       for (m=0; m<n; m++) { f1[m] = tm * g1[m];    f2[m] = tm * g2[m]; }
  1283.     else {            /* odd power ==> do one more mult */
  1284.       mult(g1,g2,f1,f2,w1,w2);    /* f1 and f2 contain the result now */
  1285.       for (m=0; m<n; m++) { f1[m] = tm * f1[m]; f2[m] = tm * f2[m]; }}
  1286.   }
  1287.   else {            /* BACKWARD FFT */
  1288.     if (l==1)            /* even power */
  1289.       for (m=0; m<n; m++) { f1[m] = g1[m]; f2[m] = g2[m]; }
  1290.     else            /* odd power ==> do one more mult */
  1291.       mult(g1,g2,f1,f2,w1,w2);    /* f1 and f2 contain the result now */
  1292.   }
  1293. }
  1294.  
  1295.  
  1296. /* CHIRP-Z-TRANSFORM (complex data)
  1297.  */
  1298.  
  1299. /* C_Array_CZT           Generalization of DFT
  1300.    ;;
  1301.    Frequency is scaled as an L/2-point DFT of the input data (zero padded to L/2).
  1302.    ;;
  1303.    phi = starting point (on unit circle) -- Range 0,1 (covers 0,2pi like DFT angle)
  1304.    rho = resolution (angular frequency spacing) -- Range 0,1 (maps 0,2pi like DFT angle)
  1305.    N = input data length
  1306.    M = output data length
  1307.    log2_L = smallest_power_of_2_ge(N+M-1)   ----
  1308.    f1,f2 contain the input data (complex).
  1309.    f1,f2,fo1,fo2,g1,g2            must be of length L
  1310.    fft_w1,fft_w2                  must be of length L/2
  1311.    czt_w1,czt_w2                  must be of length max(M,N)  ----
  1312.    ;;
  1313.    RESULT is left on f1,f2 (M complex numbers).
  1314.    */
  1315. C_Array_CZT(phi,rho, N,M,log2_L, f1,f2,fo1,fo2, g1,g2, fft_w1,fft_w2,czt_w1,czt_w2)
  1316.      double phi, rho;
  1317.      REAL *f1,*f2,*fo1,*fo2, *g1,*g2,  *fft_w1,*fft_w2,*czt_w1,*czt_w2;
  1318.      long N,M,log2_L;
  1319. { long i, maxMN, L, L2;
  1320.   void Make_CZT_Tables(), CZT_Pre_Multiply(), Make_Chirp_Filter();
  1321.   void Make_FFT_Tables(), C_Array_FFT_With_Given_Tables(), C_Array_Complex_Multiply_Into_First_One();
  1322.  
  1323.   maxMN = max(M,N);
  1324.   L = 1<<log2_L;
  1325.   L2 = L/2;
  1326.  
  1327.   CZT_Pre_Multiply(phi,rho, f1,f2, N,L);
  1328.   Make_FFT_Tables(fft_w1,fft_w2, L, 1);    /* PREPARE TABLES FOR FORWARD FFT */
  1329.   C_Array_FFT_With_Given_Tables(1, log2_L, L, f1,f2, g1,g2, fft_w1,fft_w2);
  1330.  
  1331.   Make_CZT_Tables(czt_w1,czt_w2, rho, maxMN);
  1332.   Make_Chirp_Filter(fo1,fo2, N,M,L, czt_w1,czt_w2);
  1333.   C_Array_FFT_With_Given_Tables(1, log2_L, L, fo1,fo2, g1,g2, fft_w1,fft_w2);
  1334.  
  1335.   C_Array_Complex_Multiply_Into_First_One(f1,f2, fo1,fo2, L);
  1336.   for (i=0;i<L2;i++) fft_w2[i] = (-fft_w2[i]); /* PREPARE TABLES FOR INVERSE FFT */
  1337.   C_Array_FFT_With_Given_Tables(-1, log2_L, L, f1,f2, g1,g2, fft_w1,fft_w2);
  1338.   C_Array_Complex_Multiply_Into_First_One(f1,f2, czt_w1,czt_w2, M);
  1339. }
  1340.  
  1341. void CZT_Pre_Multiply(phi,rho, f1,f2, N,L)      /* phi = starting point */
  1342.      double phi,rho; REAL *f1,*f2; long N,L;    /* this proc is two complex multiplication */
  1343. { long i;
  1344.   double tmp, A1, A2;
  1345.   rho = rho*.5;            /* To make 1/2 in exponent "(n^2)/2" */
  1346.   for (i=0;i<N;i++)
  1347.   { tmp = ((double) i);
  1348.     tmp = TWOPI * ((phi + (rho*tmp)) * tmp);
  1349.     A1 = cos(tmp);
  1350.     A2 = sin(tmp);
  1351.     tmp   =         A1*f1[i] - A2*f2[i];
  1352.     f2[i] = (REAL) (A1*f2[i] + A2*f1[i]);
  1353.     f1[i] = (REAL) tmp;
  1354.   }
  1355.   for (i=N;i<L;i++) { f1[i] = 0.0; /* zero pad */
  1356.               f2[i] = 0.0; }
  1357. }
  1358.  
  1359. void Make_Chirp_Filter(fo1,fo2, N,M,L, czt_w1,czt_w2)
  1360.      REAL *fo1,*fo2, *czt_w1,*czt_w2; long N,M,L;
  1361. { long i, L_minus_N_plus_1 = L-N+1;
  1362.   for (i=0;i<M;i++) { fo1[i] =   czt_w1[i];
  1363.               fo2[i] = - czt_w2[i]; }
  1364.   for (i=M;i<L_minus_N_plus_1;i++) { fo1[i] = 0.0; /* arbitrary region, circular convolution */
  1365.                      fo2[i] = 0.0; }
  1366.   for (i=L_minus_N_plus_1;i<L;i++) { fo1[i] =   czt_w1[(L-i)];
  1367.                      fo2[i] = - czt_w2[(L-i)]; }
  1368. }
  1369.  
  1370. void Make_CZT_Tables(czt_w1,czt_w2, rho, maxMN)              /* rho = resolution */
  1371.      double rho; REAL *czt_w1,*czt_w2; long maxMN;
  1372. { long i;
  1373.   double tmp;
  1374.   rho = rho*.5;            /* the 1/2 in the "(n^2)/2" exponent */
  1375.   for (i=0;i<maxMN;i++)
  1376.   { tmp = ((double) i);
  1377.     tmp = TWOPI * (tmp * (rho*tmp));
  1378.     czt_w1[i] = (REAL) cos(tmp);
  1379.     czt_w2[i] = (REAL) sin(tmp); }
  1380. }
  1381.  
  1382.  
  1383. long smallest_power_of_2_ge(n)
  1384.      long n;
  1385. { long i,power;
  1386.   if (n<0) { printf("\n ABORT program! smallest_pwr_of_2_ge negative argument--- %d\n", n); fflush(stdout); }
  1387.   power=0; i=1;
  1388.   while (i<n)
  1389.   { power++; i=i*2; }
  1390.   return(power);
  1391. }
  1392.  
  1393. /*  stuff not currently used
  1394.  
  1395. void CZT_Post_Multiply(f1,f2,czt_w1,czt_w2,M)
  1396.      REAL *f1,*f2,*czt_w1,*czt_w2; long M;
  1397. { long i;
  1398.   REAL tmp;
  1399.   for (i=0;i<M;i++)
  1400.   { tmp = f1[i]*czt_w1[i] - f2[i]*czt_w2[i];
  1401.     f2[i] = 2.0 * (f1[i]*czt_w2[i] + f2[i]*czt_w1[i]);
  1402.     f1[i] = 2.0 * tmp;
  1403.   }}
  1404.  
  1405. #define take_modulo_one(x, answer)  \
  1406. { long ignore_integral_part;        \
  1407.   double modf();                    \
  1408.   answer = (double) modf( ((double) x), &ignore_integral_part); }
  1409.             ^ this only works when answer is double
  1410.  
  1411. */
  1412.  
  1413.  
  1414.  
  1415. /* 2D DFT ---------------- row-column decomposition
  1416.    (3D not working yet)
  1417.  */
  1418. C_Array_2D_FFT_In_Scheme_Heap(flag, nrows, ncols, Real_Array, Imag_Array)
  1419.      long flag, nrows, ncols; REAL *Real_Array, *Imag_Array;
  1420. { long i, j;
  1421.   REAL *Temp_Array;
  1422.   REAL *f1,*f2,*g1,*g2,*w1,*w2, *Work_Here;
  1423.   long nrows_power, ncols_power, Length = nrows*ncols;
  1424.  
  1425.   if (nrows==ncols) {        /* SQUARE IMAGE, OPTIMIZE... */
  1426.     Square_Image_2D_FFT_In_Scheme_Heap(flag, nrows, Real_Array, Imag_Array);
  1427.   }
  1428.   else {            /* NOT A SQUARE IMAGE, CANNOT DO FAST_TRANSPOSE */
  1429.     /* FIRST (NCOLS-1)POINT FFTS FOR EACH ROW, THEN (NROWS-1)POINT FFTS FOR EACH COLUMN */
  1430.  
  1431.     for (ncols_power=0, i=ncols; i>1; ncols_power++) { /* FIND/CHECK POWERS OF ROWS,COLS */
  1432.       if ( (i % 2) == 1) error_bad_range_arg (2);
  1433.       i=i/2; }
  1434.     for (nrows_power=0, i=nrows; i>1; nrows_power++) {
  1435.       if ( (i % 2) == 1) error_bad_range_arg (1);
  1436.       i=i/2; }
  1437.  
  1438. #if (REAL_IS_DEFINED_DOUBLE != 0)
  1439.     ALIGN_FLOAT (Free);
  1440.     Free += 1;
  1441. #endif
  1442.     Primitive_GC_If_Needed(Length*REAL_SIZE + ((max(nrows,ncols))*3*REAL_SIZE));
  1443.     Work_Here = (REAL *) Free;
  1444.     g1 = Work_Here;
  1445.     g2 = Work_Here + ncols;
  1446.     w1 = Work_Here + (ncols<<1);
  1447.     w2 = Work_Here + (ncols<<1) + (ncols>>1);
  1448.     Make_FFT_Tables(w1,w2,ncols, flag);
  1449.     for (i=0;i<nrows;i++) {    /* ROW-WISE */
  1450.       f1 = Real_Array + (i*ncols);
  1451.       f2 = Imag_Array + (i*ncols);
  1452.       C_Array_FFT_With_Given_Tables(flag, ncols_power, ncols, f1,f2,g1,g2,w1,w2);
  1453.     }
  1454.  
  1455.     Temp_Array = Work_Here;
  1456.     Work_Here  = Temp_Array + Length;
  1457.     Image_Transpose(Real_Array, Temp_Array, nrows, ncols); /* TRANSPOSE: (1) order of frequencies. (2) read columns.*/
  1458.     Image_Transpose(Imag_Array, Real_Array, nrows, ncols);
  1459.  
  1460.     g1 = Work_Here;
  1461.     g2 = Work_Here + nrows;
  1462.     w1 = Work_Here + (nrows<<1);
  1463.     w2 = Work_Here + (nrows<<1) + (nrows>>1);
  1464.     Make_FFT_Tables(w1,w2,nrows,flag);
  1465.     for (i=0;i<ncols;i++) {    /* COLUMN-WISE */
  1466.       f1 = Temp_Array + (i*nrows); /* THIS IS REAL DATA */
  1467.       f2 = Real_Array + (i*nrows); /* THIS IS IMAG DATA */
  1468.       C_Array_FFT_With_Given_Tables(flag, nrows_power, nrows, f1,f2,g1,g2,w1,w2);
  1469.     }
  1470.  
  1471.     Image_Transpose(Real_Array, Imag_Array, ncols, nrows); /* DO FIRST THIS !!!, do not screw up Real_Data !!! */
  1472.     Image_Transpose(Temp_Array, Real_Array, ncols, nrows); /* TRANSPOSE BACK: order of frequencies. */
  1473.   }
  1474. }
  1475.  
  1476. Square_Image_2D_FFT_In_Scheme_Heap(flag, nrows, Real_Array, Imag_Array)
  1477.      long flag,nrows; REAL *Real_Array, *Imag_Array;
  1478. { REAL *f1,*f2,*g1,*g2,*w1,*w2, *Work_Here;
  1479.   long nrows_power;
  1480.   long i;
  1481.  
  1482.   for (nrows_power=0, i=nrows; i>1; nrows_power++) { /* FIND/CHECK POWERS OF ROWS */
  1483.     if ( (i % 2) == 1) error_bad_range_arg (2);
  1484.     i=i/2; }
  1485. #if (REAL_IS_DEFINED_DOUBLE != 0)
  1486.     ALIGN_FLOAT (Free);
  1487.     Free += 1;
  1488. #endif
  1489.   Primitive_GC_If_Needed(nrows*3*REAL_SIZE);
  1490.   Work_Here = (REAL *) Free;
  1491.   g1 = Work_Here;
  1492.   g2 = Work_Here + nrows;
  1493.   w1 = Work_Here + (nrows<<1);
  1494.   w2 = Work_Here + (nrows<<1) + (nrows>>1);
  1495.   Make_FFT_Tables(w1, w2, nrows, flag);    /* MAKE TABLES */
  1496.   for (i=0;i<nrows;i++) {    /* ROW-WISE */
  1497.     f1 = Real_Array + (i*nrows);
  1498.     f2 = Imag_Array + (i*nrows);
  1499.     C_Array_FFT_With_Given_Tables(flag, nrows_power, nrows, f1,f2,g1,g2,w1,w2);
  1500.   }
  1501.   Image_Fast_Transpose(Real_Array, nrows); /* MUST TRANSPOSE (1) order of frequencies. (2) read columns. */
  1502.   Image_Fast_Transpose(Imag_Array, nrows);
  1503.  
  1504.   for (i=0;i<nrows;i++) {    /* COLUMN-WISE */
  1505.     f1 = Real_Array + (i*nrows);
  1506.     f2 = Imag_Array + (i*nrows);
  1507.     C_Array_FFT_With_Given_Tables(flag, nrows_power, nrows, f1,f2,g1,g2,w1,w2);    /* ncols=nrows... Twiddles... */
  1508.   }
  1509.   Image_Fast_Transpose(Real_Array, nrows); /* TRANSPOSE BACK: order of frequencies. */
  1510.   Image_Fast_Transpose(Imag_Array, nrows);
  1511. }
  1512.  
  1513. C_Array_3D_FFT_In_Scheme_Heap(flag, ndeps, nrows, ncols, Real_Array, Imag_Array)
  1514.      long flag, ndeps, nrows, ncols; REAL *Real_Array, *Imag_Array;
  1515. { long l, m, n;
  1516.   REAL *Temp_Array;
  1517.   REAL *f1,*f2,*g1,*g2,*w1,*w2, *Work_Here;
  1518.   long ndeps_power, nrows_power, ncols_power;
  1519.  
  1520.   if ((ndeps==nrows) && (nrows==ncols)) {                                           /* CUBIC IMAGE, OPTIMIZE... */
  1521.     Cube_Space_3D_FFT_In_Scheme_Heap(flag, ndeps, Real_Array, Imag_Array);
  1522.   }
  1523.   else {
  1524.     for (ndeps_power=0, l=ndeps; l>1; ndeps_power++) {                 /* FIND/CHECK POWERS OF DEPS,ROWS,COLS */
  1525.       if ( (l % 2) == 1) error_bad_range_arg (2);
  1526.       l=l/2; }
  1527.     for (nrows_power=0, m=nrows; m>1; nrows_power++) {
  1528.       if ( (m % 2) == 1) error_bad_range_arg (2);
  1529.       m=m/2; }
  1530.     for (ncols_power=0, n=ncols; n>1; ncols_power++) {
  1531.       if ( (n % 2) == 1) error_bad_range_arg (2);
  1532.       n=n/2; }
  1533.  
  1534.     printf("3D FFT implemented only for cubic-spaces.\n");
  1535.     printf("aborted\n.");
  1536.   }
  1537. }
  1538.  
  1539. Cube_Space_3D_FFT_In_Scheme_Heap(flag, ndeps, Real_Array, Imag_Array)
  1540.      long flag, ndeps; REAL *Real_Array, *Imag_Array;
  1541. { fast long l, m, n;
  1542.   fast long ndeps_power, Surface_Length;
  1543.   fast REAL *From_Real, *From_Imag;
  1544.   fast REAL *f1,*f2,*g1,*g2,*w1,*w2, *Work_Here;
  1545.  
  1546.   for (ndeps_power=0, l=ndeps; l>1; ndeps_power++) {                 /* FIND/CHECK POWER OF NDEPS */
  1547.     if ( (l % 2) == 1) error_bad_range_arg (2);
  1548.     l=l/2; }
  1549. #if (REAL_IS_DEFINED_DOUBLE != 0)
  1550.     ALIGN_FLOAT (Free);
  1551.     Free += 1;
  1552. #endif
  1553.   Primitive_GC_If_Needed(ndeps*3*REAL_SIZE);
  1554.   Work_Here = (REAL *) Free;
  1555.   g1 = Work_Here;
  1556.   g2 = Work_Here + ndeps;
  1557.   w1 = Work_Here + (ndeps<<1);
  1558.   w2 = Work_Here + (ndeps<<1) + (ndeps>>1);
  1559.   Make_FFT_Tables(w1, w2, ndeps, flag);                      /* MAKE TABLES */
  1560.  
  1561.   Surface_Length=ndeps*ndeps;
  1562.   From_Real = Real_Array;   From_Imag = Imag_Array;
  1563.  
  1564.   for (l=0; l<ndeps; l++,From_Real+=Surface_Length,From_Imag+=Surface_Length) {       /* DEPTH-WISE */
  1565.  
  1566.     f1 = From_Real;    f2 = From_Imag;
  1567.     for (m=0; m<ndeps; m++,f1+=ndeps,f2+=ndeps) {                                     /* ROW-WISE */
  1568.       C_Array_FFT_With_Given_Tables(flag, ndeps_power, ndeps, f1,f2,g1,g2,w1,w2); }
  1569.     Image_Fast_Transpose(From_Real, ndeps);    /* MUST TRANSPOSE (1) order of frequencies. (2) read columns. */
  1570.     Image_Fast_Transpose(From_Imag, ndeps);
  1571.  
  1572.     /* ndeps=nrows=ncols, same Twiddle Tables */
  1573.  
  1574.     f1 = From_Real;    f2 = From_Imag;
  1575.     for (n=0; n<ndeps; n++,f1+=ndeps,f2+=ndeps) {                                   /* COLUMN-WISE */
  1576.       C_Array_FFT_With_Given_Tables(flag, ndeps_power, ndeps, f1,f2,g1,g2,w1,w2); }
  1577.     Image_Fast_Transpose(From_Real, ndeps);            /* TRANSPOSE BACK: order of frequencies. */
  1578.     Image_Fast_Transpose(From_Imag, ndeps);
  1579.   }
  1580. }
  1581.  
  1582.  
  1583. /*----------------------- scheme primitives ------------------------- */
  1584.  
  1585. /* Real and Imag arrays must be different.
  1586.    Arg1=1 --> forward FFT, otherwise backward.
  1587.    */
  1588.  
  1589. DEFINE_PRIMITIVE ("ARRAY-FFT!", Prim_array_fft, 3, 3, 0)
  1590. { long length, power, flag, i;
  1591.   SCHEME_OBJECT answer;
  1592.   REAL *f1,*f2,*g1,*g2,*w1,*w2;
  1593.   REAL *Work_Here;
  1594.  
  1595.   PRIMITIVE_HEADER (4);
  1596.   flag = arg_integer(1);    /* forward or backward  */
  1597.   CHECK_ARG (2, ARRAY_P);    /*      input real */
  1598.   CHECK_ARG (3, ARRAY_P);    /*      input imag */
  1599.  
  1600.   length = ARRAY_LENGTH(ARG_REF(2));
  1601.   if (length != (ARRAY_LENGTH(ARG_REF(3)))) error_bad_range_arg(2);
  1602.  
  1603.   for (power=0, i=length; i>1; power++) {
  1604.     if ( (i % 2) == 1) error_bad_range_arg(2);
  1605.     i=i/2; }
  1606.  
  1607.   f1 = ARRAY_CONTENTS(ARG_REF(2));
  1608.   f2 = ARRAY_CONTENTS(ARG_REF(3));
  1609.   if (f1==f2)  error_wrong_type_arg(2);
  1610.  
  1611. #if (REAL_IS_DEFINED_DOUBLE != 0)
  1612.     ALIGN_FLOAT (Free);
  1613.     Free += 1;
  1614. #endif
  1615.   Primitive_GC_If_Needed(length*3*REAL_SIZE);
  1616.   Work_Here = (REAL *) Free;
  1617.   g1 = Work_Here;
  1618.   g2 = Work_Here + length;
  1619.   w1 = Work_Here + (length<<1);
  1620.   w2 = Work_Here + (length<<1) + (length>>1);
  1621.  
  1622.   C_Array_FFT(flag, power, length, f1,f2,g1,g2,w1,w2);
  1623.  
  1624.   PRIMITIVE_RETURN (UNSPECIFIC);
  1625. }
  1626.  
  1627. DEFINE_PRIMITIVE ("ARRAY-CZT!", Prim_array_czt, 6,6, 0)
  1628. { double phi,rho;
  1629.   long N,M,L,  i;
  1630.   long log2_L, maxMN;
  1631.   long smallest_power_of_2_ge();
  1632.   int errcode;
  1633.   REAL *a,*b,*c,*d;
  1634.   REAL *f1,*f2,*fo1,*fo2, *g1,*g2, *fft_w1,*fft_w2,*czt_w1,*czt_w2,    *Work_Here;
  1635.  
  1636.   PRIMITIVE_HEADER (6);
  1637.   phi = (arg_real_number (1));    /* starting point [0,1]*/
  1638.   rho = (arg_real_number (2));    /* resolution [0,1] */
  1639.   CHECK_ARG (3, ARRAY_P);    /* input real */
  1640.   CHECK_ARG (4, ARRAY_P);    /* input imag */
  1641.   CHECK_ARG (5, ARRAY_P);    /* output real */
  1642.   CHECK_ARG (6, ARRAY_P);    /* output imag */
  1643.   
  1644.   a = ARRAY_CONTENTS(ARG_REF(3));
  1645.   b = ARRAY_CONTENTS(ARG_REF(4));
  1646.   c = ARRAY_CONTENTS(ARG_REF(5));
  1647.   d = ARRAY_CONTENTS(ARG_REF(6));
  1648.   
  1649.   N = ARRAY_LENGTH(ARG_REF(3));    /* N = input length */
  1650.   M = ARRAY_LENGTH(ARG_REF(5));    /* M = output length */
  1651.   if (N!=(ARRAY_LENGTH(ARG_REF(4))))    error_bad_range_arg(3);
  1652.   if (M!=(ARRAY_LENGTH(ARG_REF(6))))    error_bad_range_arg(5);
  1653.  
  1654.   if ((M+N-1) < 4)                      error_bad_range_arg(5);
  1655.   log2_L = smallest_power_of_2_ge(M+N-1);
  1656.   L  = 1<<log2_L;        /* length of intermediate computation arrays */
  1657.   maxMN =  (((M)<(N)) ? (N) : (M)); /* length of czt tables =  maximum(M,N) */
  1658.  
  1659. #if (REAL_IS_DEFINED_DOUBLE != 0)
  1660.     ALIGN_FLOAT (Free);
  1661.     Free += 1;
  1662. #endif
  1663.   Primitive_GC_If_Needed( ((7*L) + (2*maxMN)) * REAL_SIZE);
  1664.   g1  = (REAL *) Free;
  1665.   g2  = g1  + L;
  1666.   fo1 = g2  + L;
  1667.   fo2 = fo1 + L;
  1668.   fft_w1 = fo2 + L;
  1669.   fft_w2 = fft_w1 + (L/2);
  1670.   czt_w1 = fft_w2 + (L/2);
  1671.   czt_w2 = czt_w1 + maxMN;
  1672.   f1 = czt_w2 + maxMN;        /* CZT stores its results here */
  1673.   f2 = f1     + L;
  1674.  
  1675.   for (i=0; i<N; i++) { f1[i] = a[i]; /*        input data */
  1676.             f2[i] = b[i]; }
  1677.  
  1678.   C_Array_CZT(phi,rho, N,M,log2_L, f1,f2,fo1,fo2, g1,g2, fft_w1,fft_w2,czt_w1,czt_w2);
  1679.  
  1680.   for (i=0; i<M; i++) { c[i] = f1[i]; /*        results */
  1681.             d[i] = f2[i]; }
  1682.  
  1683.   PRIMITIVE_RETURN (UNSPECIFIC);
  1684. }
  1685.  
  1686. DEFINE_PRIMITIVE ("ARRAY-2D-FFT!", Prim_array_2d_fft, 5, 5, 0)
  1687. {
  1688.   PRIMITIVE_HEADER (5);
  1689.   {
  1690.     fast long nrows = (arg_integer_in_range (2, 1, 513));
  1691.     fast long ncols = (arg_integer_in_range (3, 1, 513));
  1692.     fast SCHEME_OBJECT real_image = (ARG_REF (4));
  1693.     fast SCHEME_OBJECT imag_image = (ARG_REF (5));
  1694.     CHECK_ARG (4, ARRAY_P);
  1695.     CHECK_ARG (5, ARRAY_P);
  1696.     if (real_image == imag_image)
  1697.       error_wrong_type_arg (5);
  1698.     Set_Time_Zone (Zone_Math);
  1699.   {
  1700.     long length = (ARRAY_LENGTH (real_image));
  1701.     if ((length != (ARRAY_LENGTH (imag_image))) ||
  1702.     (length != (nrows * ncols)))
  1703.       error_bad_range_arg (5);
  1704.   }
  1705.     C_Array_2D_FFT_In_Scheme_Heap
  1706.       ((arg_integer (1)),    /* flag 1=forward else backward */
  1707.        nrows,
  1708.        ncols,
  1709.        (ARRAY_CONTENTS (real_image)),
  1710.        (ARRAY_CONTENTS (imag_image)));
  1711.     PRIMITIVE_RETURN (cons (real_image, (cons (imag_image, EMPTY_LIST))));
  1712.   }
  1713. }
  1714.  
  1715. DEFINE_PRIMITIVE ("ARRAY-3D-FFT!", Prim_array_3d_fft, 6, 6, 0)
  1716. {
  1717.   PRIMITIVE_HEADER (6);
  1718.   {
  1719.     fast long ndeps = (arg_integer_in_range (2, 1, 513));
  1720.     fast long nrows = (arg_integer_in_range (3, 1, 513));
  1721.     fast long ncols = (arg_integer_in_range (4, 1, 513));
  1722.     fast SCHEME_OBJECT real_image = (ARG_REF (5));
  1723.     fast SCHEME_OBJECT imag_image = (ARG_REF (6));
  1724.     CHECK_ARG (5, ARRAY_P);
  1725.     CHECK_ARG (6, ARRAY_P);
  1726.     if (real_image == imag_image)
  1727.       error_wrong_type_arg (6);
  1728.     {
  1729.       long length = (ARRAY_LENGTH (real_image));
  1730.       if ((length != (ARRAY_LENGTH (imag_image))) ||
  1731.       (length != (ndeps * nrows * ncols)))
  1732.     error_bad_range_arg (6);
  1733.     }
  1734.     C_Array_3D_FFT_In_Scheme_Heap
  1735.       ((arg_integer (1)),
  1736.        ndeps,
  1737.        nrows,
  1738.        ncols,
  1739.        (ARRAY_CONTENTS (real_image)),
  1740.        (ARRAY_CONTENTS (imag_image)));
  1741.     PRIMITIVE_RETURN (cons (real_image, (cons (imag_image, EMPTY_LIST))));
  1742.   }
  1743. }
  1744.